LPA and SC

suppressPackageStartupMessages({
  library("here")
  library("tidyverse")
  library("MplusAutomation")
  library("gt")
  library("glue")
  library("kableExtra")
  library("misty")
  library("lavaan")
  library("AICcmodavg")
  library("nonnest2")
  library("DiagrammeR")
  library("lavaan")
  library("tidyLPA")
  library("semTools")
  library("brms")
  library("MBESS")
  library("ufs")
  library("robmed")
  library("careless")
  library("mice")
  library("psych")
  library("BayesFactor")
  library("effectsize")
  library("tidybayes")
  library("emmeans")
  library("bayesplot")
  library("patchwork")
  library("bmlm")
})

# install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))


options("max.print" = .Machine$integer.max)

# Make random things reproducible
set.seed(1234)

options(mc.cores = 4)

bayesplot_theme_set()

source(here::here("src", "functions", "funs_add_neoffi60_subscales.R"))
source(here::here("src", "functions", "funs_correct_iesr_scores.R"))
source(here::here("src", "functions", "funs_plot_job_qualification.R"))
source(here::here("src", "functions", "funs_generate_all_items_df.R"))

scale_this <- function(x) as.vector(scale(x))

sum_coding <- function(x, lvls = levels(x)) {
  # codes the first category with -1
  nlvls <- length(lvls)
  stopifnot(nlvls > 1)
  cont <- diag(nlvls)[, -nlvls, drop = FALSE]
  cont[nlvls, ] <- -1
  cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
  colnames(cont) <- lvls[-1]
  x <- factor(x, levels = lvls)
  contrasts(x) <- cont
  x
}

Get data

all_items <- generate_all_items_df()

There is a problem with IES-R, in the control group. I shift the control distribution of IES-R towards lower values.

```r
temp <- correct_iesr_scores(all_items)
all_items <- temp

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->




<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuZ2dwbG90KGFsbF9pdGVtcywgYWVzKHggPSBpZXNyX3RzLCBjb2xvdXIgPSBpc19yZXNjdWVfd29ya2VyKSkgK1xuICBnZW9tX2RlbnNpdHkoKVxuYGBgIn0= -->

```r
ggplot(all_items, aes(x = iesr_ts, colour = is_rescue_worker)) +
  geom_density()

all_items |> 
  ggplot(aes(y=iesr_ts, fill=is_rescue_worker)) + 
    geom_boxplot() +
    scale_y_log10()

all_items$y <- all_items$iesr_ts + 0.001

all_items |> 
  group_by(is_rescue_worker) |> 
  summarize(
    # gm_iesr = geometric.mean(y, na.rm = TRUE),
    avg_iesr = mean(y, na.rm = TRUE)
  )

There are some cases in which the RWs of Toscana and Lombardia did not have the proper qualification. They are assigned to the category of RWs.

all_items$idx <- 1:nrow(all_items)
all_items$is_rescue_worker <- as.character(all_items$is_rescue_worker)
all_items$is_rescue_worker <- ifelse(
  all_items$idx < 746, "Si", all_items$is_rescue_worker
)
all_items$idx <- NULL

all_items$commeetee_location <- ifelse(
  all_items$is_rescue_worker == "No", "None", all_items$red_cross_commeetee_location
) |> 
  factor()

all_items$commeetee_location <- ifelse(
  all_items$is_rescue_worker == "No" & all_items$age < 25, "students", all_items$commeetee_location
)

all_items$commeetee <- ifelse(
  all_items$is_rescue_worker == "No" & all_items$age >= 25, "community_sample",
  all_items$commeetee_location
)

all_items$commeetee <- 
  ifelse(is.na(all_items$commeetee), "community_sample", all_items$commeetee) |> 
  factor()

sum(is.na(all_items$commeetee))

[1] 0

```r
```r
```r
# all_items <- all_items %>%
#   mutate(job_qualification = case_when(
#     is_rescue_worker==\Si\ & job_qualification == \non_rescue_worker\ ~ \team_member\,
#     TRUE ~ job_qualification)) 

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



## IES-R as a function of group

iesr_ts | trunc(lb = 0) ~ is_rescue_worker + (1 | commeetee),


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTAgPC0gYnJtKFxuICBiZihpZXNyX3RzIH4gaXNfcmVzY3VlX3dvcmtlciksXG4gIGZhbWlseSA9IGh1cmRsZV9nYW1tYSgpLFxuICBkYXRhID0gYWxsX2l0ZW1zLFxuICBiYWNrZW5kID0gXCJjbWRzdGFuclwiXG4gICMgYWxnb3JpdGhtID0gXCJtZWFuZmllbGRcIlxuKVxuYGBgIn0= -->

```r
m0 <- brm(
  bf(iesr_ts ~ is_rescue_worker),
  family = hurdle_gamma(),
  data = all_items,
  backend = "cmdstanr"
  # algorithm = "meanfield"
)

pp_check(m0)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m0)

 Family: hurdle_gamma 
  Links: mu = log; shape = identity; hu = identity 
Formula: iesr_ts ~ is_rescue_worker 
   Data: all_items (Number of observations: 1068) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              2.69      0.05     2.60     2.79 1.00     4431     2990
is_rescue_workerSi     0.26      0.06     0.15     0.37 1.00     4401     3236

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     1.41      0.06     1.30     1.52 1.00     4553     2820
hu        0.05      0.01     0.04     0.06 1.00     4772     3170

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m0, "is_rescue_worker"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  all_items$ies_ts[all_items$is_rescue_worker == "Si"], 
  all_items$ies_ts[all_items$is_rescue_worker == "No"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.28      | [0.15, 0.41]

- Estimated using pooled SD.

Supported families are: ‘acat’, ‘asym_laplace’, ‘bernoulli’, ‘beta’, ‘beta_binomial’, ‘binomial’, ‘categorical’, ‘com_poisson’, ‘cox’, ‘cratio’, ‘cumulative’, ‘custom’, ‘dirichlet’, ‘dirichlet2’, ‘discrete_weibull’, ‘exgaussian’, ‘exponential’, ‘frechet’, ‘gamma’, ‘gaussian’, ‘gen_extreme_value’, ‘geometric’, ‘hurdle_cumulative’, ‘hurdle_gamma’, ‘hurdle_lognormal’, ‘hurdle_negbinomial’, ‘hurdle_poisson’, ‘info’, ‘inverse.gaussian’, ‘logistic_normal’, ‘lognormal’, ‘multinomial’, ‘negbinomial’, ‘negbinomial2’, ‘poisson’, ‘shifted_lognormal’, ‘skew_normal’, ‘sratio’, ‘student’, ‘von_mises’, ‘weibull’, ‘wiener’, ‘zero_inflated_asym_laplace’, ‘zero_inflated_beta’, ‘zero_inflated_beta_binomial’, ‘zero_inflated_binomial’, ‘zero_inflated_negbinomial’, ‘zero_inflated_poisson’, ‘zero_one_inflated_beta’

The sk, ch, mi sub-scales are coded so that high values indicate high self-compassion levels. The sj, is, oi sub-scales are coded so that high values indicate low self-compassion levels.

The ts_sc score has been computed by reversing the coding of the items of the sj, is, oi sub-scales (so that they indicate the absence of self-judgment, absence of isolation, absence of over-identification).

scs_subscales <- with(all_items, data.frame(sk, ch, mi, sj, is, oi, scs_ts))
cor(scs_subscales) |> round(2)

          sk    ch    mi    sj    is    oi scs_ts
sk      1.00  0.52  0.58 -0.39 -0.28 -0.24   0.71
ch      0.52  1.00  0.49 -0.01 -0.03 -0.04   0.45
mi      0.58  0.49  1.00 -0.19 -0.33 -0.35   0.66
sj     -0.39 -0.01 -0.19  1.00  0.67  0.66  -0.75
is     -0.28 -0.03 -0.33  0.67  1.00  0.80  -0.78
oi     -0.24 -0.04 -0.35  0.66  0.80  1.00  -0.78
scs_ts  0.71  0.45  0.66 -0.75 -0.78 -0.78   1.00

COPE scale

In the COPE scale only two factors are identified.

all_items$pos_reinterpretation <- with(all_items, cope_1 + cope_29 + cope_38 + cope_59)
all_items$mental_disengagement <- with(all_items, cope_2 + cope_16 + cope_31 + cope_43) 
all_items$venting <- with(all_items, cope_3 + cope_17 + cope_28 + cope_46) 
all_items$seeking_instrumental_support <- with(all_items, cope_4 + cope_14 + cope_30 + cope_45) 
all_items$active_coping <- with(all_items, cope_5 + cope_25 + cope_47 + cope_58)  
all_items$denial <- with(all_items, cope_6 + cope_27 + cope_40 + cope_57) 
all_items$religion <- with(all_items, cope_7 + cope_18 + cope_48 + cope_60) 
all_items$humor <- with(all_items, cope_8 + cope_20 + cope_36 + cope_50) 
all_items$behavioral_disengagement <- with(all_items, cope_9 + cope_24 + cope_37 + cope_51) 
all_items$restraint <- with(all_items, cope_10 + cope_22 + cope_41 + cope_49) 
all_items$seeking_emotional_support <- with(all_items, cope_11 + cope_23 + cope_34 + cope_52) 
all_items$substance_use <- with(all_items, cope_12 + cope_26 + cope_35 + cope_53) 
all_items$acceptance <- with(all_items, cope_13 + cope_21 + cope_44 + cope_54) 
all_items$suppr_competing_activities <- with(all_items, cope_15 + cope_33 + cope_42 + cope_55) 
all_items$planning <- with(all_items, cope_19 + cope_32 + cope_39 + cope_56) 

Create COPE sub-scales scores using all items – note that SEM analyses suggest to drop some of the items.

all_items$active_coping <- with(
  all_items, pos_reinterpretation + active_coping +
  suppr_competing_activities + planning + restraint + 
    seeking_instrumental_support + acceptance
)

all_items$avoidance_coping <- with(
  all_items, mental_disengagement + denial + humor +
  behavioral_disengagement + substance_use + religion 
)

all_items$soc_emo_coping <- with(
  all_items, seeking_instrumental_support +
  seeking_emotional_support + venting
)

Self-compassion scale

plot(density(all_items$scs_ts))

fit_1 <- brm(
  bf(
    scs_ts ~ is_rescue_worker,
    sigma ~ is_rescue_worker
    ),
  family = student(),
  backend = "cmdstanr",
  data = all_items
)

pp_check(fit_1)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

me <- conditional_effects(
  fit_1, "is_rescue_worker"
)
plot(me, points = FALSE)

summary(fit_1)

 Family: student 
  Links: mu = identity; sigma = log; nu = identity 
Formula: scs_ts ~ is_rescue_worker 
         sigma ~ is_rescue_worker
   Data: all_items (Number of observations: 1068) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                   74.11      1.02    72.15    76.10 1.00     3862     2674
sigma_Intercept              2.88      0.04     2.80     2.97 1.00     4140     2724
is_rescue_workerSi           7.49      1.19     5.17     9.78 1.00     4050     2963
sigma_is_rescue_workerSi    -0.13      0.05    -0.23    -0.03 1.00     4027     2960

Family Specific Parameters: 
   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu    38.40     15.75    16.97    77.94 1.00     4228     3115

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

BFt <- BayesFactor::ttestBF(
  all_items$scs_ts[all_items$is_rescue_worker == "Si"], 
  all_items$scs_ts[all_items$is_rescue_worker == "No"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.42      | [0.29, 0.56]

- Estimated using pooled SD.

rw_df <- all_items |> 
  dplyr::filter(is_rescue_worker == "Si")

rw_df <- rw_df %>%
  mutate(job_qualification = case_when(
    job_qualification == "non_rescue_worker" ~ "team_member",
    TRUE ~ job_qualification))

LPA

lpa_scales <- c(
  "is_rescue_worker",
  "neuroticism", "extraversion", "openness", "agreeableness", "conscientiousness",
  "active_coping", "avoidance_coping", "soc_emo_coping",
  "iesr_ts",
  # "avoiding", "intrusivity", "hyperarousal",
  # "sk", "ch", "mi", "sj", "is", "oi",
  # "pos_sc",
  # "neg_sc",
  # "ts_sc",
  "mpss_tot"
  # "ptgi_total_score"
  # "relating_to_others",
  # "new_possibilities",
  # "personal_strength",
  # "appreciation_of_life",
  # "spirituality"
)

lpa_rw_df <- subset(rw_df, select=lpa_scales) 

# lpa_rw_df <- subset(all_items, select=lpa_scales) |> 
#   dplyr::filter(is_rescue_worker == "Si")

lpa_rw_df$is_rescue_worker <- NULL

head(lpa_rw_df)

lpa_rw_df %>% 
  scale() %>%
  estimate_profiles(1:6,
    variances = c("equal", "varying"),
    covariances = c("zero", "varying")
    #package = "MplusAutomation"
  )

tidyLPA analysis using mclust: 

 Model Classes AIC      BIC      Entropy prob_min prob_max n_min n_max BLRT_p
 1     1       21342.45 21434.88 1.00    1.00     1.00     1.00  1.00        
 1     2       20740.50 20883.76 0.71    0.87     0.94     0.36  0.64  0.01  
 1     3       20459.15 20653.25 0.75    0.82     0.93     0.18  0.57  0.01  
 1     4       20314.46 20559.39 0.72    0.84     0.84     0.11  0.41  0.01  
 1     5       20232.20 20527.97 0.71    0.78     0.84     0.09  0.32  0.01  
 1     6       20161.53 20508.13 0.74    0.72     0.89     0.05  0.32  0.01  
 6     1       19882.54 20182.93 1.00    1.00     1.00     1.00  1.00        
 6     2       19509.14 20114.55 0.70    0.89     0.94     0.44  0.56  0.01  
 6     3       19415.46 20325.88 0.73    0.87     0.90     0.26  0.38  0.01  
 6     4       19396.37 20611.80 0.72    0.82     0.89     0.23  0.29  0.19  
 6     5       19359.41 20879.85 0.77    0.84     0.90     0.10  0.28  0.07  
 6     6       19321.00 21146.46 0.81    0.81     0.99     0.08  0.23  0.04  

lpa_rw_df %>% 
  scale() %>%
  estimate_profiles(1:6,
    variances = c("equal", "varying"),
    covariances = c("zero", "varying")
    # package = "MplusAutomation"
  ) %>% 
  compare_solutions(statistics = c("AIC", "BIC"))

m2 <- lpa_rw_df %>%
  scale() %>%
  estimate_profiles(2,
    variances = "varying",
    covariances = "varying",
    package = "MplusAutomation"
  )

m2_plot <- lpa_rw_df %>%
  scale() %>%
  estimate_profiles(2,
    variances = "varying",
    covariances = "varying",
    package = "MplusAutomation"
    ) %>%
    plot_profiles(add_line = TRUE, rawdata= FALSE, bw = FALSE)

Profile 2: dysfunctional Profile 1: adaptive

get_estimates(m2)

out <- get_data(m2)
lpa_rw_df$lpa_class <- out$Class

table(
  lpa_rw_df$lpa_class
)

  1   2 
409 342 

table(
  lpa_rw_df$lpa_class, rw_df$job_qualification
)

   
    driver team_leader team_member
  1     96         161         152
  2     59         138         145

lpa_rw_df$class <- factor(lpa_rw_df$lpa_class)
summary(lpa_rw_df$class)

  1   2 
409 342 

rw_df$class <- lpa_rw_df$class
rw_df$profile <- lpa_rw_df$class
rw_df$profile <- ifelse(
  rw_df$profile == "2", "maladaptive", "adaptive"
)
rw_df$profile <- factor(rw_df$profile)
# Reorder the levels of the 'profile' factor
rw_df$profile <- factor(rw_df$profile, levels = c("maladaptive", "adaptive"))

m1 <- brm(
  bf(scs_ts ~ profile),
  family = skew_normal(),
  data = rw_df,
  init = 0.1,
  backend = "cmdstanr",
  adapt_delta = 0.9
)

pp_check(m1)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

summary(m1)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: scs_ts ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          74.55      0.83    72.90    76.14 1.00     3443     2621
profileadaptive    12.64      1.15    10.39    14.92 1.00     3493     2699

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    14.90      0.38    14.19    15.68 1.00     4123     2947
alpha    -0.70      0.59    -1.53     0.58 1.00     2187     3498

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m1, "profile"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  rw_df$scs_ts[rw_df$class == 1], 
  rw_df$scs_ts[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt, type = "d")

Cohen's d |       95% CI
------------------------
0.87      | [0.71, 1.02]

- Estimated using pooled SD.

m1 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  labs(x = "LPA Class", y = "SCS Score", title = "Rescue Workers") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 83, label = "Bayesian Cohen's d = 0.89\n 95% CI [0.73, 1.04]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

names(all_items)

Self-judgment

m2 <- brm(
  bf(sj ~ profile),
  data = rw_df,
  family = student,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m2)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m2 <- loo(m2)
plot(loo_m2)

summary(m2)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: sj ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          17.03      0.23    16.58    17.48 1.00     3615     2904
profileadaptive    -3.61      0.31    -4.20    -2.99 1.00     3754     2741

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.27      0.12     4.04     4.52 1.00     2939     2589
nu       46.35     16.93    21.54    86.78 1.00     3388     2924

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m2, "profile"
)
plot(me, points = FALSE)

emmeans(m2, specs = pairwise ~ profile)

$emmeans
 profile     emmean lower.HPD upper.HPD
 maladaptive   17.0      16.6      17.5
 adaptive      13.4      13.0      13.8

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast               estimate lower.HPD upper.HPD
 maladaptive - adaptive     3.61      3.01      4.21

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$sj[rw_df$profile == "adaptive"], 
  rw_df$sj[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |         95% CI
--------------------------
-0.81     | [-0.97, -0.66]

- Estimated using pooled SD.

p2 <- m2 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Self-Judgment", title = "A") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.82\n 95% CI [0.67, 0.97]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p2

Isolation

m3 <- brm(
  bf(is ~ profile),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.95
)

pp_check(m3)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m3 <- loo(m3)
plot(loo_m3)

summary(m3)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: is ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          12.29      0.25    11.78    12.80 1.00     2470     2490
profileadaptive    -3.38      0.34    -4.05    -2.69 1.00     2358     2099

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.06      0.12     3.84     4.32 1.00     2494     2593
alpha     2.31      0.65     1.13     3.65 1.00     1915     2113

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m3, "profile"
)
plot(me, points = FALSE)

emmeans(m3, specs = pairwise ~ profile)

$emmeans
 profile     emmean lower.HPD upper.HPD
 maladaptive  12.29     11.80     12.81
 adaptive      8.91      8.49      9.28

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast               estimate lower.HPD upper.HPD
 maladaptive - adaptive     3.38      2.74      4.09

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$is[rw_df$profile == "adaptive"], 
  rw_df$is[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |         95% CI
--------------------------
-0.93     | [-1.08, -0.78]

- Estimated using pooled SD.

p3 <- m3 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Isolation", title = "B") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 10.2, label = "Bayesian Cohen's d = 0.94\n 95% CI [0.79, 1.10]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p3

Over-identification

m4 <- brm(
  bf(oi ~ profile),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m4)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m4 <- loo(m4)
plot(loo_m4)

summary(m4)

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: oi ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          10.96      0.24    10.48    11.42 1.00     1492     1549
profileadaptive    -2.51      0.34    -3.16    -1.85 1.00     1366     1537

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.67      0.12     3.44     3.92 1.00     1369     1758
alpha     3.59      0.95     2.16     5.73 1.00     1253     1261

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m4, "profile"
)
plot(me, points = FALSE)

emmeans(m4, specs = pairwise ~ profile)

$emmeans
 profile     emmean lower.HPD upper.HPD
 maladaptive  10.97     10.48     11.42
 adaptive      8.44      8.04      8.79

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast               estimate lower.HPD upper.HPD
 maladaptive - adaptive     2.52      1.86      3.16

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$oi[rw_df$profile == "adaptive"], 
  rw_df$oi[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |         95% CI
--------------------------
-0.96     | [-1.11, -0.81]

- Estimated using pooled SD.

p4 <- m4 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Over-Identification", title = "C") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 9.4, label = "Bayesian Cohen's d = 0.97\n 95% CI [0.82, 1.12]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p4

Self-kindness

m5 <- brm(
  bf(sk ~ profile),
  family = student(),
  data = rw_df,
  backend = "cmdstanr"
)

pp_check(m5)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m5 <- loo(m5)
plot(loo_m5)

summary(m5)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: sk ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          13.05      0.23    12.60    13.52 1.00     3991     3028
profileadaptive     1.25      0.31     0.63     1.86 1.00     4251     2936

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.21      0.12     3.98     4.45 1.00     3510     2614
nu       41.31     16.98    17.09    83.01 1.00     3245     2884

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m5, "profile"
)
plot(me, points = FALSE)

emmeans(m5, specs = pairwise ~ profile)

$emmeans
 profile     emmean lower.HPD upper.HPD
 maladaptive   13.1      12.6      13.5
 adaptive      14.3      13.9      14.7

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast               estimate lower.HPD upper.HPD
 maladaptive - adaptive    -1.26     -1.85    -0.625

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$sk[rw_df$profile == "adaptive"], 
  rw_df$sk[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.28      | [0.13, 0.42]

- Estimated using pooled SD.

p5 <- m5 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Self-Kindness", title = "D") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.28\n 95% CI [0.14, 0.43]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p5

Common humanity

m6 <- brm(
  bf(ch ~ profile),
  family = student(),
  data = rw_df,
  backend = "cmdstanr"
)

pp_check(m6)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m6 <- loo(m6)
plot(loo_m6)

summary(m6)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: ch ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          11.59      0.18    11.23    11.96 1.00     4291     3143
profileadaptive     0.01      0.24    -0.47     0.48 1.00     4425     3056

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.28      0.09     3.10     3.46 1.00     4625     2797
nu       43.03     16.67    19.04    83.45 1.00     4370     3134

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

me <- conditional_effects(
  m6, "profile"
)
plot(me, points = FALSE)

BFt <- BayesFactor::ttestBF(
  rw_df$ch[rw_df$profile == "adaptive"], 
  rw_df$ch[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |        95% CI
-------------------------
-5.07e-04 | [-0.14, 0.14]

- Estimated using pooled SD.

p6 <- m6 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Common-Humanity", title = "E") +
  papaja::theme_apa() + 
  annotate("text", x = 1.5, y = 12.2, label = "Bayesian Cohen's d = 0.00\n 95% CI [-0.14, 0.14]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p6

Mindfulness

m7 <- brm(
  bf(mi ~ profile),
  family = student(),
  data = rw_df,
  backend = "cmdstanr"
)

Compiling Stan program...

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

-

\

|

/

 

Start sampling
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: student_t_lpdf: Scale parameter is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1 
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: student_t_lpdf: Scale parameter is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1 
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2 
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: student_t_lpdf: Scale parameter is 0, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2 
Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 3 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 3 
Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 3 Exception: student_t_lpdf: Scale parameter is 0, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 75, column 4 to column 48)
Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 3 
Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 4 Exception: gamma_lpdf: Random variable is inf, but must be positive finite! (in '/var/folders/cy/4xdvhqx966nggmk95hsnyzc40000gn/T/RtmpFKK0iZ/model-9ac975dd8c9f.stan', line 67, column 2 to column 66)
Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 4 

pp_check(m7)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

loo_m7 <- loo(m7)
plot(loo_m7)

summary(m7)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: mi ~ profile 
   Data: rw_df (Number of observations: 751) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept          12.74      0.17    12.41    13.08 1.01     3843     2838
profileadaptive     1.01      0.23     0.58     1.46 1.00     3804     2789

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.03      0.09     2.86     3.21 1.00     3965     2848
nu       36.26     15.13    15.29    72.29 1.00     3682     3429

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

emmeans(m7, specs = pairwise ~ profile)

$emmeans
 profile     emmean lower.HPD upper.HPD
 maladaptive   12.7      12.4      13.1
 adaptive      13.8      13.5      14.1

Point estimate displayed: median 
HPD interval probability: 0.95 

$contrasts
 contrast               estimate lower.HPD upper.HPD
 maladaptive - adaptive    -1.01     -1.44    -0.557

Point estimate displayed: median 
HPD interval probability: 0.95 

BFt <- BayesFactor::ttestBF(
  rw_df$mi[rw_df$profile == "adaptive"], 
  rw_df$mi[rw_df$profile == "maladaptive"],
  paired = FALSE
)
effectsize(BFt)

Cohen's d |       95% CI
------------------------
0.32      | [0.18, 0.46]

- Estimated using pooled SD.

me <- conditional_effects(
  m7, "profile"
)
plot(me, points = FALSE)

p7 <- m7 %>%
  emmeans( ~ profile) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = profile, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "SCS Mindfulness", title = "F") +
  papaja::theme_apa() + 
  annotate("text", x = 1, y = 13.8, label = "Bayesian Cohen's d = 0.32\n 95% CI [0.17, 0.46]")

Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").

p7

fig_scs <- (p2 | p3 | p4) /
(p5 | p6 | p7)

out <- fig_scs + plot_annotation(
  title = 'SCS Subscales as a Function of LPA Profile'
  # subtitle = 'Rescue Workers group'
  # caption = 'Disclaimer: None of these plots are insightful'
)
ggsave("scs_subscales_lpa.pdf", width = 35, height = 20, units = "cm")

print(out)

IES-R

m10 <- brm(
  bf(ies_ts ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m10)

summary(m10)

p10 <- m10 %>%
  emmeans( ~ class) %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = class, y = .value)) +
  geom_eye() +
  stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
  # facet_grid(~ wool) +
  # theme_light()
  scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
  labs(x = "LPA Class", y = "Impact of Event Scale - Revised (IES-R)") +
  # papaja::theme_apa() + 
  annotate("text", x = 1, y = 17, label = "Bayesian Cohen's d = 1.34\n 95% CI [1.18, 1.50]")
p10

BFt <- BayesFactor::ttestBF(
  rw_df$ies_ts[rw_df$class == 1], 
  rw_df$ies_ts[rw_df$class == 2],
  paired = FALSE
)
effectsize(BFt)

emmeans(m10 , specs = pairwise ~ class)

m11 <- brm(
  bf(ptgi_total_score | trunc(lb = 0) ~ class),
  family = student(),
  data = rw_df,
  backend = "cmdstanr"
)

pp_check(m11)

summary(m11)

Job qualification

rw_df$job_qualification <- ifelse(
  rw_df$job_qualification == "non_rescue_worker", "team_member", 
  rw_df$job_qualification
) 

m12 <- brm(
  bf(ies_ts ~ job_qualification),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m12)

summary(m12)

me <- conditional_effects(
  m12, "job_qualification"
)
plot(me, points = FALSE)

emmeans(m12 , specs = pairwise ~ job_qualification)

m13 <- brm(
  bf(scs_ts ~ job_qualification * class),
  family = gaussian(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m13)

me <- conditional_effects(
  m13, "job_qualification:class"
)
plot(me, points = FALSE)

summary(m13)

emmeans(m13 , specs = pairwise ~ job_qualification*class)

m14 <- brm(
  bf(ptgi_total_score ~ ies_total_score | (1 | comme\)),
  family = gaussian(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

summary(m14)

bayes_R2(m14)

conditional_effects(m14, "is:class")

SEM

mydf <- data.frame(
  scs = scale(rw_df$scs_ts),
  class = ifelse(rw_df$class == 1, 0.0, 1.0),
  ptgi = scale(rw_df$ptgi_total_score),
  psc = scale(rw_df$sk + rw_df$ch + rw_df$mi),
  nsc = scale(rw_df$sj + rw_df$oi + rw_df$is),
  commettee = rw_df$red_cross_commeetee_location,
  id = 1:nrow(rw_df)
)

mydf <- mydf[complete.cases(mydf), ]

# Impute NAs with the mode.
# mydf$commettee[is.na(mydf$commettee)] <- "Comitato di Groane"
# summary(factor(mydf$commettee))

f1 <- bf(scs ~ class, family = skew_normal())
f2 <- bf(ptgi ~ scs + class, family = skew_normal())
mod <- brm(
  f1 + f2 + set_rescor(FALSE), 
  data = mydf, 
  cores = 4, 
  refresh = 0,
  init = 0.1,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

bayestestR::mediation(mod, mediator = "scs", ci = 0.95, method = "SPI")

pp_check(mod, resp = "ptgi")

pp_check(mod, resp = "scs")

summary(mod)

fit <- mlm(d = mydf, 
           id = "commettee",
           x = "class",
           m = "scs",
           y = "ptgi",
           iter = 2000, 
           cores = 4)

mlm_path_plot(fit, level = .95, text = T,
              xlab = "Resilience\nProfile",
              mlab = "Self\nCompassion",
              ylab = "PTG", digits = 2)

Causal Mediation Analysis for Stan Model

Treatment: class Mediator : scs Response : ptgi

Effect | Estimate | 95% SPI

Direct Effect (ADE) | -0.300 | [-0.464, -0.140] Indirect Effect (ACME) | 0.121 | [ 0.053, 0.187] Mediator Effect | 0.154 | [ 0.069, 0.230] Total Effect | -0.178 | [-0.325, -0.027]

m16 <- brm(
  bf(ptgi_total_score ~ class),
  family = skew_normal(),
  data = rw_df,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m16)

summary(m16)

bayes_R2(m16)

m17 <- brm(
  bf(ptgi_total_score ~ is_rescue_worker),
  family = skew_normal(),
  data = all_items,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

summary(m17)

rw_df$rate_of_activity_num <- as.integer(rw_df$rate_of_activity_num)

temp <- rw_df[1:746, ]

temp$activity <- cut(
  temp$rate_of_activity_num,
  breaks = c(-1, 0, 1, 2),    
  labels = c("Low", "Medium", "High"),   
  include.lowest = TRUE,    
  ordered_result = TRUE
)   

m19 <- brm(
  bf(activity ~ class),
  family=cumulative("logit"),
  data = temp,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

pp_check(m19)

summary(m19)

m20 <- brm(
  bf(class ~ last_training_num * scs_ts),
  family=bernoulli(),
  data = temp,
  backend = "cmdstanr",
  adapt_delta = 0.99
)

summary(m20)

conditional_effects(m20, "scs_ts:last_training_num")

table(rw_df$class, rw_df$last_training_num)